GLM Part6

Preparations

Load the necessary libraries

library(car)       #for regression diagnostics
library(broom)     #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot)    #for outputs
library(knitr)     #for kable
library(effects)   #for partial effects plots
library(emmeans)   #for estimating marginal means
library(ggeffects) #for plotting marginal means
library(MASS)      #for glm.nb
library(MuMIn)     #for AICc
library(tidyverse) #for data wrangling
library(modelr)    #for auxillary modelling functions
library(DHARMa)    #for residual diagnostics plots
library(performance) #for residuals diagnostics
library(see)         #for plotting residuals
library(patchwork) #grid of plots
library(scales)    #for more scales

Scenario

An ecologist studying a rocky shore at Phillip Island, in southeastern Australia, was interested in how clumps of intertidal mussels are maintained (Quinn 1988). In particular, he wanted to know how densities of adult mussels affected recruitment of young individuals from the plankton. As with most marine invertebrates, recruitment is highly patchy in time, so he expected to find seasonal variation, and the interaction between season and density - whether effects of adult mussel density vary across seasons - was the aspect of most interest.

The data were collected from four seasons, and with two densities of adult mussels. The experiment consisted of clumps of adult mussels attached to the rocks. These clumps were then brought back to the laboratory, and the number of baby mussels recorded. There were 3-6 replicate clumps for each density and season combination.

Format of quinn.csv data files

SEASON	DENSITY	RECRUITS	SQRTRECRUITS	GROUP
Spring	Low	15	3.87	SpringLow
..	..	..	..	..
Spring	High	11	3.32	SpringHigh
..	..	..	..	..
Summer	Low	21	4.58	SummerLow
..	..	..	..	..
Summer	High	34	5.83	SummerHigh
..	..	..	..	..
Autumn	Low	14	3.74	AutumnLow
..	..	..	..	..

SEASON	Categorical listing of Season in which mussel clumps were collected ­ independent variable
DENSITY	Categorical listing of the density of mussels within mussel clump ­ independent variable
RECRUITS	The number of mussel recruits ­ response variable
SQRTRECRUITS	Square root transformation of RECRUITS - needed to meet the test assumptions
GROUPS	Categorical listing of Season/Density combinations - used for checking ANOVA assumptions

Mussel

Read in the data

quinn = read_csv('../public/data/quinn.csv', trim_ws=TRUE)
glimpse(quinn)

## Rows: 42
## Columns: 5
## $ SEASON       <chr> "Spring", "Spring", "Spring", "Spring", "Spring", "Sprin…
## $ DENSITY      <chr> "Low", "Low", "Low", "Low", "Low", "High", "High", "High…
## $ RECRUITS     <dbl> 15, 10, 13, 13, 5, 11, 10, 15, 10, 13, 1, 21, 31, 21, 18…
## $ SQRTRECRUITS <dbl> 3.872983, 3.162278, 3.605551, 3.605551, 2.236068, 3.3166…
## $ GROUP        <chr> "SpringLow", "SpringLow", "SpringLow", "SpringLow", "Spr…

summary(quinn)

##     SEASON            DENSITY             RECRUITS      SQRTRECRUITS  
##  Length:42          Length:42          Min.   : 0.00   Min.   :0.000  
##  Class :character   Class :character   1st Qu.: 9.25   1st Qu.:3.041  
##  Mode  :character   Mode  :character   Median :13.50   Median :3.674  
##                                        Mean   :18.33   Mean   :3.871  
##                                        3rd Qu.:21.75   3rd Qu.:4.663  
##                                        Max.   :69.00   Max.   :8.307  
##     GROUP          
##  Length:42         
##  Class :character  
##  Mode  :character  
##                    
##                    
## 

Since we intend to model both SEASON and DENSITY as categorical variables, we need to explicitly declair them as factors.

quinn = quinn %>% mutate(SEASON = factor(SEASON, levels=c('Spring', 'Summer', 'Autumn', 'Winter')),
                         DENSITY = factor(DENSITY))

Exploratory data analysis

Model formula: \[ \begin{align} y_i &\sim{} \mathcal{Pois}(\lambda_i)\\ ln(\mu_i) &= \boldsymbol{\beta} \bf{X_i}\\[1em] \end{align} \]

where \(\boldsymbol{\beta}\) is a vector of effects parameters and \(\bf{X}\) is a model matrix representing the intercept and effects of season, density and their interaction on mussel recruitment.

ggplot(quinn, aes(y=RECRUITS, x=SEASON, fill=DENSITY)) +
     geom_boxplot()

Conclusions:

• there is clear evidence of non-homogeneity of variance
• specifically, there is evidence that the variance is related to the mean in that boxplots that are lower on the y-axis (low mean) also have lower variance (shorter boxplots)
• this might be expected for count data and we might consider that a Poisson distribution (which assumes that mean and variance are equal - and thus related in a very specific way).

Lets mimic the effect of using a log link, by using log scaled y-axis.

ggplot(quinn, aes(y=RECRUITS, x=SEASON, fill=DENSITY)) +
  geom_boxplot() +
  scale_y_log10()

Conclusions:

• that is an improvement

Fit the model

We actually have a couple of initial options:

1. log transform the response variable (RECRUITS) and fit against a Gaussian distribution. If we do so, we would need to ammend the response for cases when the response is zero (as the log of 0 is not defined). We could simply add 1 to each count before log transformed.

2. fit against a Poisson distribution with a log link

log-transformed and Gaussian

I will include this, just for an illustration on how to fit such a model, we will not persue it further.

quinn.glmG <- glm(log(RECRUITS+1) ~ DENSITY*SEASON, data=quinn, family=gaussian)

Poisson

quinn.glm <- glm(RECRUITS ~ DENSITY*SEASON, data=quinn,
                  family=poisson(link='log'))

Model validation

autoplot

autoplot(quinn.glm,which=1:6)

Conclusions:

• most of these diagnostics seem ok, however, there is one observation (#34) that has a very high residual (and thus Cook’s D)

Performance model checking

performance::check_model(quinn.glm)

performance::check_overdispersion(quinn.glm)

## # Overdispersion test
## 
##        dispersion ratio =   3.309
##   Pearson's Chi-Squared = 112.497
##                 p-value = < 0.001

performance::check_zeroinflation(quinn.glm)

## # Check for zero-inflation
## 
##    Observed zeros: 2
##   Predicted zeros: 0
##             Ratio: 0.00

Conclusions:

• there is evidence of overdispersion

DHARMa residuals

quinn.resid = simulateResiduals(quinn.glm,  plot=TRUE)

testResiduals(quinn.resid)

## $uniformity
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.2008, p-value = 0.06761
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.7813, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 6, observations = 42, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.04761905
## sample estimates:
## outlier frequency (expected: 0.00880952380952381 ) 
##                                          0.1428571

## $uniformity
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.2008, p-value = 0.06761
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.7813, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 6, observations = 42, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.04761905
## sample estimates:
## outlier frequency (expected: 0.00880952380952381 ) 
##                                          0.1428571

testDispersion(quinn.resid)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.7813, p-value < 2.2e-16
## alternative hypothesis: two.sided

testZeroInflation(quinn.resid)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 10.87, p-value = 0.016
## alternative hypothesis: two.sided

#testTemporalAutocorrelation(quinn.glm1)

Conclusions:

• there is evidence of overdispersion and outliers
• there is also some evidence of zero-inflation

Simple dispersion

## goodness of fit
1-pchisq(quinn.glm$deviance, df=quinn.glm$df.residual)

## [1] 2.457479e-12

## any evidence of overdispersion
quinn.glm$deviance/quinn.glm$df.residual

## [1] 3.677366

Conclusions:

• there is evidence of over-dispersion

Simple zero-inflation

quinn %>% group_by(SEASON, DENSITY) %>%
  summarise(Zeros=sum(RECRUITS==0), 
            Prop=Zeros/n(),
            Mean=mean(RECRUITS))

x=rpois(100000,lambda=2.67)
tab.1 = table(x==0)
tab.1/sum(tab.1)

## 
##   FALSE    TRUE 
## 0.93188 0.06812

## OR,  over the entire data
## is this due to excessive zeros (zero-inflation)
tab=table(quinn$RECRUITS==0)
tab/sum(tab)

## 
##      FALSE       TRUE 
## 0.95238095 0.04761905

## 5% is not many.. so it cant be zero-inflated
## how many 0's would we expect from a poisson distribution with a mean similar to our mean
mean(quinn$RECRUITS)

## [1] 18.33333

x=rpois(100000,lambda=mean(quinn$RECRUITS))
tab.1 = table(x==0)
tab.1/sum(tab.1)

## 
## FALSE 
##     1

Conclusion:

• although at the level of the entire dataset, there is no evidence for excessive zeros, if we explore this separately within each SEASON by DENSITY combination, we see that for the Winter/Low density group, the proportion of zeros is very high.

Different model

In light of the discovery that the Poisson model was overdispersed, we will explore a negative binomial model. Rather than assume that the variance is equal to the mean (dispersion of 1), the dispersion parameter is estimated. Negative binomial models are also able to account for some level of zero-inflation.

Fit model

library(MASS)
quinn.nb = glm.nb(RECRUITS ~ DENSITY*SEASON, data=quinn)

Model validation

autoplot

autoplot(quinn.nb,which=1:6)

Performance model checking

performance::check_model(quinn.nb)

DHARMa residuals

quinn.resid = simulateResiduals(quinn.nb,  plot=TRUE)

testResiduals(quinn.resid)

## $uniformity
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.11213, p-value = 0.626
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.89834, p-value = 0.624
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 42, p-value = 0.66
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.04761905
## sample estimates:
## outlier frequency (expected: 0.0102380952380952 ) 
##                                        0.02380952

## $uniformity
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.11213, p-value = 0.626
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.89834, p-value = 0.624
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 42, p-value = 0.66
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.04761905
## sample estimates:
## outlier frequency (expected: 0.0102380952380952 ) 
##                                        0.02380952

testDispersion(quinn.resid)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.89834, p-value = 0.624
## alternative hypothesis: two.sided

testZeroInflation(quinn.resid)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 6.1728, p-value = 0.104
## alternative hypothesis: two.sided

Simple dispersion

## goodness of fit
1-pchisq(quinn.nb$deviance, df=quinn.nb$df.residual)

## [1] 0.01311451

## any evidence of overdispersion
quinn.nb$deviance/quinn.nb$df.residual

## [1] 1.614218

Model comparison

AICc(quinn.glm, quinn.nb)

Conclusions:

• on the basis of AICc, the negative binomial model would be considered better than the Poisson model.

Conclusions:

• there is no evidence that the negative binomial model is overdispersed or zero-inflated
• all other diagnostics are also improved

Parial plots

plot_model

plot_model(quinn.nb, type='eff',  terms=c('SEASON', 'DENSITY'))

allEffects

plot(allEffects(quinn.nb),multiline=TRUE, ci.style='bar')

plot(allEffects(quinn.nb),multiline=TRUE, ci.style='bar', type='link')

ggpredict

ggpredict(quinn.nb, c('SEASON', 'DENSITY')) %>% plot()

ggemmeans

ggemmeans(quinn.nb,  ~SEASON*DENSITY) %>% plot()

Model investigation / hypothesis testing

summary

Lets start with the estimated coeffcicients on the link (log) scale

summary(quinn.nb)

## 
## Call:
## glm.nb(formula = RECRUITS ~ DENSITY * SEASON, data = quinn, init.theta = 9.022467857, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8704  -0.8274   0.0000   0.4999   2.1866  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               2.3026     0.1875  12.283  < 2e-16 ***
## DENSITYLow                0.1133     0.2742   0.413  0.67934    
## SEASONSummer              1.5721     0.2389   6.581 4.69e-11 ***
## SEASONAutumn              0.6763     0.2492   2.714  0.00664 ** 
## SEASONWinter             -0.5680     0.2881  -1.971  0.04870 *  
## DENSITYLow:SEASONSummer  -0.8970     0.3509  -2.556  0.01059 *  
## DENSITYLow:SEASONAutumn  -0.1881     0.3788  -0.496  0.61955    
## DENSITYLow:SEASONWinter  -0.8671     0.5338  -1.624  0.10432    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(9.0225) family taken to be 1)
## 
##     Null deviance: 183.269  on 41  degrees of freedom
## Residual deviance:  54.883  on 34  degrees of freedom
## AIC: 293.09
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  9.02 
##           Std. Err.:  3.69 
## 
##  2 x log-likelihood:  -275.095

Conclusions:

• the intercept represents the estimated mean of the first combination of Season (Spring) and Density (High). On the link scale this is 2.3
• the difference between Low and High adult density in spring is 0.11, although this is not significant
• the difference between Spring and Summer for High adult density is 1.57
• the difference between Spring and Autumn for High adult density is 0.68
• the difference between Spring and Winter for High adult density is -0.57
• if there were no interactions, we would expect the Low density Summer recruitment to be the additive of the main effects (Low and Summer). However, the modelled mean is 0.9 less than the additive effects would have expected. This value is significantly different to 0, indicating that there is evidence that the density effect in Summer is different to that in Spring.
• the density effect in Autum and Winter were not found to be significantly different from what you would expect from an additive model.

tidy

tidy(quinn.nb, conf.int=TRUE)

Now if we exponentiate to put these estimates on the response scale:

tidy(quinn.nb, conf.int=TRUE, exponentiate=TRUE)

Conclusions:

• the intercept represents the estimated mean of the first combination of Season (Spring) and Density (High). On the response scale this is 10
• there is a 0.11 fold difference between High and Low adult density in spring, although this is not significant
• there is a 1.57 fold difference between Summer and Spring for High adult density
• there is a 0.68 fold difference between Autumn and Spring for High adult density
• there is a -0.57 fold difference between Winter and Spring for High adult density
  
• the modelled mean for Summer Low density is 0.9 fold different (less than half) that we would have expected in the absence of interactions
• the density effect in Autum and Winter were not found to be significantly different from what you would expect from an additive model.

If we compare the above to the overdispersed poisson, we see that the estimates are the same, but that the standard errors are underestimated and hence the confidence intervals are narrower (for overdispersed model)

tidy(quinn.glm, conf.int=TRUE, exponentiate=TRUE)

Predictions

In order to tease appart the significant interaction(s), it might be useful to explore the effect of Density separately within each Season.

quinn.nb %>% emmeans(pairwise~DENSITY|SEASON)

## $emmeans
## SEASON = Spring:
##  DENSITY emmean    SE  df asymp.LCL asymp.UCL
##  High     2.303 0.187 Inf     1.935      2.67
##  Low      2.416 0.200 Inf     2.024      2.81
## 
## SEASON = Summer:
##  DENSITY emmean    SE  df asymp.LCL asymp.UCL
##  High     3.875 0.148 Inf     3.584      4.16
##  Low      3.091 0.161 Inf     2.775      3.41
## 
## SEASON = Autumn:
##  DENSITY emmean    SE  df asymp.LCL asymp.UCL
##  High     2.979 0.164 Inf     2.657      3.30
##  Low      2.904 0.203 Inf     2.505      3.30
## 
## SEASON = Winter:
##  DENSITY emmean    SE  df asymp.LCL asymp.UCL
##  High     1.735 0.219 Inf     1.306      2.16
##  Low      0.981 0.402 Inf     0.192      1.77
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
## SEASON = Spring:
##  contrast   estimate    SE  df z.ratio p.value
##  High - Low  -0.1133 0.274 Inf -0.413  0.6793 
## 
## SEASON = Summer:
##  contrast   estimate    SE  df z.ratio p.value
##  High - Low   0.7836 0.219 Inf  3.577  0.0003 
## 
## SEASON = Autumn:
##  contrast   estimate    SE  df z.ratio p.value
##  High - Low   0.0748 0.261 Inf  0.286  0.7749 
## 
## SEASON = Winter:
##  contrast   estimate    SE  df z.ratio p.value
##  High - Low   0.7538 0.458 Inf  1.646  0.0999 
## 
## Results are given on the log (not the response) scale.

Conclusions:

• there is no evidence of an effect of Density in either Spring, Autumn or Winter
• in Summer, recruitment is 0.78 higher (on a log scale) in High Density populations than Low Density populations.

Or we could express these on a response scale.

quinn.nb %>% emmeans(pairwise~DENSITY|SEASON, type='response')

## $emmeans
## SEASON = Spring:
##  DENSITY response   SE  df asymp.LCL asymp.UCL
##  High       10.00 1.87 Inf      6.93     14.44
##  Low        11.20 2.24 Inf      7.57     16.58
## 
## SEASON = Summer:
##  DENSITY response   SE  df asymp.LCL asymp.UCL
##  High       48.17 7.13 Inf     36.03     64.39
##  Low        22.00 3.55 Inf     16.03     30.19
## 
## SEASON = Autumn:
##  DENSITY response   SE  df asymp.LCL asymp.UCL
##  High       19.67 3.23 Inf     14.26     27.13
##  Low        18.25 3.71 Inf     12.25     27.19
## 
## SEASON = Winter:
##  DENSITY response   SE  df asymp.LCL asymp.UCL
##  High        5.67 1.24 Inf      3.69      8.70
##  Low         2.67 1.07 Inf      1.21      5.87
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
## SEASON = Spring:
##  contrast   ratio    SE  df z.ratio p.value
##  High / Low 0.893 0.245 Inf -0.413  0.6793 
## 
## SEASON = Summer:
##  contrast   ratio    SE  df z.ratio p.value
##  High / Low 2.189 0.480 Inf  3.577  0.0003 
## 
## SEASON = Autumn:
##  contrast   ratio    SE  df z.ratio p.value
##  High / Low 1.078 0.282 Inf  0.286  0.7749 
## 
## SEASON = Winter:
##  contrast   ratio    SE  df z.ratio p.value
##  High / Low 2.125 0.973 Inf  1.646  0.0999 
## 
## Tests are performed on the log scale

Conclusions:

• there is no evidence of an effect of Density in either Spring, Autumn or Winter
• in Summer, recruitment is 2.19 fold higher in High Density populations than Low Density populations.

quinn.nb %>% emmeans(pairwise~DENSITY|SEASON,  type='response') %>%
  confint

## $emmeans
## SEASON = Spring:
##  DENSITY response   SE  df asymp.LCL asymp.UCL
##  High       10.00 1.87 Inf      6.93     14.44
##  Low        11.20 2.24 Inf      7.57     16.58
## 
## SEASON = Summer:
##  DENSITY response   SE  df asymp.LCL asymp.UCL
##  High       48.17 7.13 Inf     36.03     64.39
##  Low        22.00 3.55 Inf     16.03     30.19
## 
## SEASON = Autumn:
##  DENSITY response   SE  df asymp.LCL asymp.UCL
##  High       19.67 3.23 Inf     14.26     27.13
##  Low        18.25 3.71 Inf     12.25     27.19
## 
## SEASON = Winter:
##  DENSITY response   SE  df asymp.LCL asymp.UCL
##  High        5.67 1.24 Inf      3.69      8.70
##  Low         2.67 1.07 Inf      1.21      5.87
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
## SEASON = Spring:
##  contrast   ratio    SE  df asymp.LCL asymp.UCL
##  High / Low 0.893 0.245 Inf     0.522      1.53
## 
## SEASON = Summer:
##  contrast   ratio    SE  df asymp.LCL asymp.UCL
##  High / Low 2.189 0.480 Inf     1.425      3.36
## 
## SEASON = Autumn:
##  contrast   ratio    SE  df asymp.LCL asymp.UCL
##  High / Low 1.078 0.282 Inf     0.646      1.80
## 
## SEASON = Winter:
##  contrast   ratio    SE  df asymp.LCL asymp.UCL
##  High / Low 2.125 0.973 Inf     0.866      5.22
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

It might be useful to capture these pairwise contrasts so that we can graph them as a caterpillar plot.

In the following, I will present the effects on a log (based 2) axis. Such a scale presents doubling and halving (etc) equidistant from 1.

eff <- (quinn.nb %>%
        emmeans(pairwise~DENSITY|SEASON,  type='response') %>%
        confint)$contrasts %>%
               as.data.frame

eff %>%
  ggplot(aes(y=ratio, x=SEASON)) +
  geom_pointrange(aes(ymin=asymp.LCL, ymax=asymp.UCL)) +
  geom_hline(yintercept=1, linetype='dashed') +
  scale_x_discrete(name='') +
  scale_y_continuous(name='Density effect (High vs Low)', trans=scales::log2_trans(),
                     breaks=scales::breaks_log(base=2)) +
  coord_flip(ylim=c(0.25, 4)) +
  theme_classic()

Summary figures

newdata = emmeans(quinn.nb, ~DENSITY|SEASON, type='response') %>% as.data.frame
head(newdata)

ggplot(newdata, aes(y=response, x=SEASON, fill=DENSITY)) +
  geom_pointrange(aes(ymin=asymp.LCL, ymax=asymp.UCL, shape=DENSITY),
                  position=position_dodge(width=0.2)) +
  theme_classic() +
  theme(axis.title.x=element_blank(),
        legend.position=c(0.01,1), legend.justification = c(0,1)) +
  annotate(geom='text', x='Summer', y=70, label='*', size=7) +
  scale_shape_manual(values=c(21, 22))

library(pscl)
quinn.zip <- zeroinfl(RECRUITS ~ DENSITY*SEASON | 1, data=quinn,  dist='poisson')
#quinn.resid <- simulateResiduals(quinn.zip,  plot=TRUE)
summary(quinn.zip)

## 
## Call:
## zeroinfl(formula = RECRUITS ~ DENSITY * SEASON | 1, data = quinn, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4510 -0.8645  0.1262  0.7225  2.8711 
## 
## Count model coefficients (poisson with log link):
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               2.3026     0.1291  17.836  < 2e-16 ***
## DENSITYLow                0.1133     0.1858   0.610  0.54191    
## SEASONSummer              1.5721     0.1419  11.081  < 2e-16 ***
## SEASONAutumn              0.6763     0.1586   4.266 1.99e-05 ***
## SEASONWinter             -0.5680     0.2147  -2.646  0.00815 ** 
## DENSITYLow:SEASONSummer  -0.8970     0.2134  -4.202 2.64e-05 ***
## DENSITYLow:SEASONAutumn  -0.1881     0.2381  -0.790  0.42957    
## DENSITYLow:SEASONWinter   0.2165     0.4534   0.478  0.63300    
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.0037     0.7305  -4.112 3.92e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 13 
## Log-likelihood: -151.3 on 9 Df

#tidy(quinn.zip,  conf.int=TRUE, exponentiate = TRUE)
exp(-3.0037)

## [1] 0.0496032

quinn.zip1 <- zeroinfl(RECRUITS ~ DENSITY*SEASON | SEASON, data=quinn,  dist='poisson')
                                        #quinn.resid <- simulateResiduals(quinn.zip)
summary(quinn.zip1)

## 
## Call:
## zeroinfl(formula = RECRUITS ~ DENSITY * SEASON | SEASON, data = quinn, 
##     dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5327 -1.0591  0.1537  0.9216  3.6831 
## 
## Count model coefficients (poisson with log link):
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               2.3026     0.1291  17.836  < 2e-16 ***
## DENSITYLow                0.1133     0.1858   0.610  0.54191    
## SEASONSummer              1.5721     0.1419  11.081  < 2e-16 ***
## SEASONAutumn              0.6763     0.1586   4.266 1.99e-05 ***
## SEASONWinter             -0.5680     0.2147  -2.646  0.00815 ** 
## DENSITYLow:SEASONSummer  -0.8970     0.2134  -4.202 2.64e-05 ***
## DENSITYLow:SEASONAutumn  -0.1881     0.2381  -0.790  0.42958    
## DENSITYLow:SEASONWinter   0.2291     0.4375   0.524  0.60045    
## 
## Zero-inflation model coefficients (binomial with logit link):
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -2.157e+01  1.454e+04  -0.001    0.999
## SEASONSummer -2.543e-08  2.013e+04   0.000    1.000
## SEASONAutumn -2.119e-08  2.107e+04   0.000    1.000
## SEASONWinter  2.031e+01  1.454e+04   0.001    0.999
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 13 
## Log-likelihood:  -148 on 12 Df

exp(-3.0037)

## [1] 0.0496032

quinn.zinb <- zeroinfl(RECRUITS ~ DENSITY*SEASON | 1, data=quinn,  dist='negbin')
AICc(quinn.zip,  quinn.zinb)

summary(quinn.zinb)

## 
## Call:
## zeroinfl(formula = RECRUITS ~ DENSITY * SEASON | 1, data = quinn, dist = "negbin")
## 
## Pearson residuals:
##        Min         1Q     Median         3Q        Max 
## -1.981e+00 -7.511e-01 -1.522e-05  5.289e-01  2.869e+00 
## 
## Count model coefficients (negbin with log link):
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               2.3026     0.1875  12.283  < 2e-16 ***
## DENSITYLow                0.1133     0.2742   0.413  0.67933    
## SEASONSummer              1.5721     0.2389   6.580 4.69e-11 ***
## SEASONAutumn              0.6763     0.2492   2.714  0.00664 ** 
## SEASONWinter             -0.5680     0.2881  -1.971  0.04869 *  
## DENSITYLow:SEASONSummer  -0.8969     0.3509  -2.556  0.01059 *  
## DENSITYLow:SEASONAutumn  -0.1881     0.3788  -0.496  0.61957    
## DENSITYLow:SEASONWinter  -0.8671     0.5338  -1.624  0.10432    
## Log(theta)                2.1997     0.4091   5.378 7.55e-08 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -15.29     453.38  -0.034    0.973
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Theta = 9.0223 
## Number of iterations in BFGS optimization: 43 
## Log-likelihood: -137.5 on 10 Df

exp(-15.29)

## [1] 2.288956e-07

References

Quinn, G. P. 1988. “Ecology of the Intertidal Pulmonate Limpet Siphonaria Diemenensis Quoy et Gaimard. II Reproductive Patterns and Energetics.” Journalofexperimentalmarinebiologyandecology 117: 137–56.